      program generate_RTF
C
C     This program must be compiled with (eg. GNU/Linux):
C
C        g77 -o genrtf  genrtf.f
C     or
C        fort77 -o genrtf  genrtf.f
C
C NB: Please keep the above commands at the top, so we can use
C     head genrtf.f command
C
C======================================================
C     Usage:
C            genrtf [-p] [-a] [-Q] [-q <int>] [-c] [-x] [-m] [-P] [-19] file
C
C====================================================== 
C     Program that
C     transforms PDB (and other formats) file into CHARMM input.
C     It creates .CRD file and also makes RTF for the residues in
C     the system, which is not in standard CHARMM distrubutions. It
C     also takes care of non AA molecules and assigns simple parameters
C     so HBUILD can be performed on the system.
C
C     HISTORY:
C
C
C     1/24/00   Version 1.0 
C               Initial version made from previous anim.f program
C               Which was build for analysis of Gaussian(TM) minimizations
C
C     12/9/00   Version 1.5
C             - many new features added
C             - Command line parameters parser added
C             - Cleanup of the code - make as much in subroutines
C                                     as possible
C             - Add XYZ and CRD format
C
C     2/4/04  - Add MOL format
C
C     4/10/01  Version 2.0
C             - First public release
C
C     7/14/04  Version 3.0
C             - Supports more input formats
C
C
C     X/YY/01  Version M.M
C             - Connect CHARMM to molden
C
C==================================================================
C
C     TODO:
C           There are lot of TODO labels in the code below - somebody fix it!
C           This is placement for general TODO list:
C
C           All the CHARMM tricks related to the generation of RTF,PSF,PARAM
C           are imlemented in version 2.0. But there is still a lot of room
C           for improvements:
C
C           - Add reading of GXX Z-matrix, include Toncka's program
C
C           - Write also .atm file for light. Facilitates to make nice movies
C             from MOLDEN trajectories (reads frames in XYZ format)
C
C           - After reading PDB file format, do the CRD file format. If HBUILD
C             cleans excessive hydrogen, we already have a final solution for
C             hydrogen addition ????
C
C           - Bumming the program. Eg on GNU/Linux see: jargon bum
C                                      (jargon is hackers' dictionary)
C
C           - Make automatic recognition of atom valences - so addition of
C             hydrogens is not so stupid any more
C
C           - Make automagic parameter assignement (can be better then MMFF
C             or connect to MMFF ???)
C
C           - Assignement of charges from QM calculations
C
C           - Run this program through the whole PDB database.
C             and watch for 'SELPRN>' how many atoms didn't make it.
C             Make the report in this header, so others can see the
C             performance, an possibly fix it! Write also a paper for
C             this. Maybe Journal of Comp. Chem??? (Better Dusanka's journal!)
C             Also do it in parallel
C             se we can make some ads for VRANA project.
C             Paper ideas:
C                - calculate contact surfaces and do the cluster analysis!!!!
C     
C
C
C==================================================================
C
C     CHECK:
C            There are some CHECK labels in the code below. These are
C            ususally in the places we are not sure if the program
C            can deal with all the situations in a correct manner...
C
C==================================================================
C     Originally designed and coded by Milan Hodoscek at
C     National Institute of Chemistry, Ljubljana, Slovenia and
C     National Institutes of Health in 1990
C
C     Code maintainer:      Milan Hodoscek, NIC Slovenia, NIH USA
C     Next maintainer:      <somebody-that-takes-over>
C                           
C     Further development:  Gregor Mlinsek, NIC Slovenia
C                           H. Lee Woodcock III, NIH, USA
C                           Benjamin (Tim) Miller, NIH, USA
C
C     COPYRIGHT: None for now
C     LICENSE:   GPL for standalone version
C                If ever build into CHARMM then follow CHARMM rules
C
C     If build into charmm the command name could be grtf. Also if
C     conflict with GenRTF arises then make this the program name
C     too.
C==================================================================

      implicit none
C     
C     Description of variables:
C==================================================================
C name(s)  max    actual Description
C          size   size   
C                        
C maxat      1      1    Maximum number of atoms (parameter)
C natom      1      1    Number of atoms in the system
C x,y,z    maxat  natom  Coordinates
C name     maxat  natom  Atom names
C resn     maxat  natom  Residue name for each atom
C chain    maxat  natom  Chain/segment name for each atom
C resid    maxat  natom  ID of residue for each atom
C ichain     1      1    Number of chain/segments in the system
C nchain   maxat  ichain Pointers to atom numbers separating the chains
C iseg       1      1    Number of chain/segments in the system
C nseg     maxat  iseg   Names of segments - same as ichain
C ires       1      1    Number of residues
C seq      maxat  ires   Names of residues
C type     maxat  natom  Name of the type: XXXX for peptide ATOMs,
C                                          First character of HETATOMs,
C                                          First 2 characters for IONS... (C*4)
C nmass      1      1    Number of new (HETATOM) types
C mass     maxat  nmass  Masses of new (HETATOM) types
C resnum     1      1    Number of residues (the entire system)
C resinx   maxat  resnum Pointers to atom numbers in each residue
C lseq     maxat  iseg   Pointers to residue numbers in each segment 
C maxpat     1      1    Maximum number of bond patches (parameter)
C ipat       1      1    Number of bond patches
C pcomm    maxpat ipat   Commands for inter residue bond patches (C*32)
C topfile    1      1    Name of topology file:from CHARMM distribution(C*255)
C parfile    1      1    Name of parameter file:from CHARMM distribution(C*255)
C ltopf      1      1    Number of characters in topfile
C lparf      1      1    Number of characters in parfile
C pdb        1      1    First 4 character from the input coordiante file (C*4)
C cresid     1      1    temp. variable for assignement of residue limits (C*4)
C cresn      1      1    temp. variable for assignement of residue limits (C*4)
C cseg       1      1    temp. variable for assignement of chain limits (C*4)
C argc       1      1    Number of command line parameters
C argv       6    argc   Command line parameters (C*80)
C iname      1      1    Resulting CHARMM input script file (C*80)
C xname      1      1    Resulting CHARMM coordinate file (CRD file) (C*80)
C tname      1      1    trajectory name - not used in this program (C*80)
C yname      1      1    name of XYZ file - not used in this program (C*80)
C tmp        1      1    Temporary string (C*80)
C chafil     1      1    Temporary buffer string for BOND commands (C*80)
C numbon     1      1    Number of bonds on HETATOM (for adding hydrogen)
C nadd       1      1    Number of added hydrogens on a given HETATOM
C parm       1      1    Logical variable - do we build parameter file template
C namgen     1      1    Counter to generate the names of atoms (H) to be added
C i,j,k      1      1    Temporary loop indexes
C nfil       1      1    Number of file parameters
C filnam     1      1    file name from command line parameters (C*100)
C comlin     1      1    command line parameters (C*100) 
C pdbflg     1      1    flag, if we want to unscramble PDB file
C charmm     1      1    flag to read CHARMM .crd file
C maxtit     1      1    Maximum number of title lines
C title      1   maxtit  Storage for CHARMM .crd titles (C*132)
C ntitle     1      1    How many lines are in the title
C resiad maxresad iresad names for HETATM residues (C*4)
C
C Tim Miller: 6/30/2008 allow segname renaming
C rnmeseg    ?      ?    Shall we rename the segment?
C segname    ?      ?    Name of segment to rename to (C*6)
C
C     For 20,000 atoms only !!!
      integer maxat
      parameter (maxat=20000)
C
      real(8) zero, one, two
      parameter (zero=0.0d0,one=1.0d0,two=2.0d0)
C
      real*8 x(maxat), y(maxat), z(maxat), wmain(maxat)
      character*4 name(maxat),resn(maxat),chain(maxat),resid(maxat)
      character*4 pdb,cresn,cresid,seq(maxat),cseg,nseg(maxat)
      character*4 newch(maxat), symb(maxat)
      logical rnmeseg
      character*6 segname
C
      character*100 comlin,filnam
      character*80 iname, tname, xname, yname,tmp,chafil
      integer argc,nfil,basenm
C
      integer maxtit,ntitle
      parameter (maxtit=60)
      character*132 title(maxtit)
      integer ltopf,lparf
      character*255 topfile,parfile
      integer maxpat
      parameter (maxpat=1000)
      character*32 pcomm(maxpat)
      integer maxresad,iresad
      parameter (maxresad=200)
      character*4 resiad(maxresad)
      logical resinew,qnew,newflag
      integer ipat,numbon,nadd,lcc
      character*2 cc
C      
      logical parm, qxyz, q19, qmol, spdb, qcon
      integer i, jlen, natom,nmass,nato,ix
      integer resnum,resinx(maxat),namgen
      integer j, k, l, ires, ion, lseq(maxat), iseg
      character*4 type(maxat), mass(maxat)
      real*8 r1, r2, requ, dist, ww
      integer ibflag, indx, ibres, ieres, ichain, nchain(maxat)
      integer ibond(maxat),jbond(maxat),nbonds
C
      integer maxyzres
      parameter (maxyzres=20000)
C
C     These are for command line parameters:
      integer qmflag,icharg,charmm,qcflag
      integer addhfl,pdbflg,xyzgrp
C
C     These are functions:
      integer islen,iargc
      real*8 dij, setmas, setrad
C
C
C     TODO:
C     Make this optional and get the info about it from outside
C     of the scope of this program:
C        getenv() or command line parameters
C
C     Also fill the array of known residues from RTF file. Now we have
C     fixed list in topres() function.
      topfile='/usr/local/charmming/toppar/top_all22_prot.inp'
      parfile='/usr/local/charmming/toppar/par_all22_prot.inp'
      ltopf=islen(topfile)
      lparf=islen(parfile)
      rnmeseg=.false.
C
      argc = iargc()
      if (argc.lt.1) then
         write(*,*)'Usage: ',
     $   'genrtf [-p] [-a] [-Q] [-C] [-q <int>] [-c] [-x] [-m] [-19] ',
     $           '[-s <segid>] file'
         write(*,*)'   -p   add template for parameters'
         write(*,*)'   -a   add hydrogens'
         write(*,*)'   -Q   add QM/MM with Gamess-US'
         write(*,*)'   -C   add QM/MM with Q-Chem'
         write(*,*)'   -q   total charge for QM'
         write(*,*)'   -c   input is CHARMM CRD'
         write(*,*)'   -x   input is XYZ'
         write(*,*)'   -m   input is MOL'
         write(*,*)'   -P   input is simple PDB'
         write(*,*)'   -19  use param19'
         write(*,*)'   -s   rename the resulting segment ID'
         write(*,*)'   -v   display the version'
         stop
      end if
C
CC HPUX      call igetarg(1,argv(1),80)
CC The rest: call getarg(1,argv(1))
C
C     Do we include all parameters needed in the system
C     as a template to build new ones
      parm=.false.
      q19=.false.
C
C     QM stuff:
      qmflag=0
      icharg=1
      qcflag=0
C
C     Do we need to add hydrogen
      addhfl=0
C
C     Is the file CHARMM coordinate file?
      charmm = 0
      pdbflg=1
C
C     for simple coordinate file formats: XYZ, MOL, simple PDB...
      qxyz = .false.
      xyzgrp=0
      qmol = .false.
      spdb=.false.
C     if there is a connection table in PDB?
      qcon=.false.
C
      nfil=0
      i = 0
 111  i=i+1
      if (i.le.argc) then
         call getarg(i,comlin)
         if (comlin(1:1).eq.'-') then
C
C           This is a place for command line options
C====
            if(comlin(2:2).eq.'v') then
               write(*,*)' GenRTF version 3.3' 
               stop
            endif 

            if(comlin(2:2).eq.'p') then
               parm=.true.
               write(*,*)'Parameters will be included'
            endif
C----
            if(comlin(2:2).eq.'x') then
               qxyz=.true.
               charmm=0
               pdbflg=0
            endif
C----
            if(comlin(2:2).eq.'m') then
               qmol=.true.
               charmm=0
               pdbflg=0
            endif
C----
            if(comlin(2:2).eq.'P') then
               spdb=.true.
               charmm=0
               pdbflg=0
            endif
C----
            if(comlin(2:2).eq.'c') then
               charmm=1
               pdbflg=0
            endif
C----
            if(comlin(2:2).eq.'Q') qmflag=1
C----
            if(comlin(2:2).eq.'C') qcflag=1
C----
            if(comlin(2:2).eq.'s') then
               i=i+1
               if (i.gt.argc) then
                  write(*,*)'Charge parameter needs a value!'
                  stop 'Wrong command'
               endif
               call getarg(i,comlin)
               read(comlin,*)segname
               rnmeseg=.true.
               write(*,*)'Renaming segment to ', segname 
            endif
C----
            if(comlin(2:2).eq.'q') then
               i=i+1
               if (i.gt.argc) then
                  write(*,*)'Charge parameter needs a value!'
                  stop 'Wrong command'
               endif
               call getarg(i,comlin)
               read(comlin,*)icharg
               write(*,*)'Charge on the QM system is: ',icharg
               goto 111
            endif
C----
            if(comlin(2:2).eq.'a')addhfl=1
C----
            if(comlin(2:2).eq.'1')q19=.true.
C----
         else
            nfil=nfil+1
            filnam(1:100)=comlin(1:100)
         endif
         goto 111
      endif
C====
      if(nfil.ne.1) then
         write(*,*)'Too many or no file parameters!'
         stop
      endif
C
C     We use only one file argument and skip file extension
C     and make .i and .t files
C
C     TODO: The following I/O block needs some cleaning
C
      basenm = index(filnam,'.')
      if (basenm.gt.0) then
         write(*,*)basenm,'dotted'
         iname=filnam(1:basenm)//'inp'
         tname=filnam(1:basenm)//'t'
         xname=filnam(1:basenm)//'crd'
         yname=filnam(1:basenm)//'xyz'
         basenm=basenm-1
         pdb=filnam(1:min(basenm,4))
      else
         jlen=islen(filnam)
         write(*,*)basenm,jlen,'dotless'
         write(*,'(1x,20i4)')(filnam(j:j),j=1,80)
         iname=filnam(1:jlen)//'.inp'
         tname=filnam(1:jlen)//'.t'
         xname=filnam(1:jlen)//'.crd'
         yname=filnam(1:jlen)//'.xyz'
         pdb=filnam(1:min(4,jlen))
         basenm=jlen
      end if
C

      open(unit=1, file=filnam, status='old')
      open(unit=2, file=iname, status='unknown')
      open(unit=3, file=tname, status='unknown')
      open(unit=4, file=xname, status='unknown')
c      open(unit=14, file=yname, status='unknown')
C
      write(*,*)'input = ',filnam
      write(*,*)'PDB = ',pdb
C
C
C     Now read the PDB file
C     Ignore everything except: ATOM, HETATOM,TER lines!
C
C
      call clrstr(maxat,chain)
      natom = 0
      ion=1
      ichain=1
      nchain(ichain)=1
C
C     READ seaction
C     We can read several formats
C
C     The default is reading PDB file
C
      if(pdbflg.eq.1) then
         call readpdb(1,15,ichain,natom,maxat,
     $        name,resn,chain,resid,x,y,z,wmain,newch,nchain,
     $        filnam(1:basenm),type,q19)
      endif
C
C     Read CHARMM .crd file: 
C
      if(charmm.eq.1) then
         call readcrd(1,15,14,filnam(1:basenm),ichain,
     $     x,y,z,wmain,title,ntitle,natom,maxat,
     $     resn,name,chain,resid,type,nchain)
      endif
C
C     Read XYZ file:
C
      if(qxyz) then
         call readxyz(1,4,15,natom,maxat,filnam(1:basenm),
     $        x,y,z,wmain,symb,name,resn,resid,chain,type,maxyzres)
      endif
C
C     Read MOL file:
C
      if(qmol) then
         call readmol(1,4,15,natom,nbonds,maxat,filnam(1:basenm),
     $        x,y,z,wmain,symb,ibond,jbond,name,resn,resid,chain,
     $        type,maxyzres)
      endif
C
C     Read simple PDB file + check for CONECT:
C
      if(spdb) then
         call readspdb(1,4,15,natom,nbonds,maxat,filnam(1:basenm),
     $        x,y,z,wmain,symb,ibond,jbond,name,resn,resid,chain,
     $        type,maxyzres,qcon)
      endif
C
C     End of main READ section
C
C     Change the types for real IONS:
C
      chain(natom+1)='FIN '
      nchain(ichain+1)=natom
C
      ion=1
      do i = 1, ichain
         if ((nchain(i+1)-nchain(i)).eq.1) then
            call settyp(2,name(nchain(i)),type(nchain(i)))
            if(ion.lt.10) then
               write(chain(nchain(i)),'(''ION'',i1)')ion
            elseif (ion.lt.100) then
               write(chain(nchain(i)),'(''IO'',i2)')ion
            endif
            ion=ion+1
         endif
      enddo
C
      write(*,*)'Number of atoms: ',natom
      if(qmol)write(*,*)'Number of bonds: ',nbonds
      write(*,*)'Number of chains/segments =',ichain
      write(*,*)'nchain(i)=',(nchain(i),i=1,ichain+1)
C
C     Fix for termini
C
      if(.not.q19)call terfix(natom,name,resid,type,chain)
C         
C     Build RTF
C============================
C
C     Fill the mass array with distinct types
C     and count them in nmass
C
C     Make MASS commnads
C
      write(2,'(''* '',a,
     $     '': automagically generated from coordinate file'')')pdb
      write(2,'(a)')'*'
      write(2,'(/''bomlev -5''/)')
      write(2,'(a,a)')'open read unit 1 card name ',topfile(1:ltopf)
      write(2,'(a)')'read rtf card unit 1'
      write(2,*)
      write(2,'(a)')'read rtf card append'
      write(2,'(a)')'* RTF'
      write(2,'(a)')'*'
      write(2,'(a)')'   28'
      write(2,*)
C
      !nmass = 1
      !mass(nmass)='H   '
      do i=1,natom
        !write(*,*)type(i)
        if (type(i).ne.'XXXX') then
           if (nmass.eq.0) then
              nmass=nmass+1
              mass(nmass)=type(i)
           else
              do j=1, nmass
                 if(type(i).eq.mass(j)) goto 888
              enddo
              nmass=nmass+1
              mass(nmass)=type(i)
 888          continue
           endif
        endif
      enddo
C
      do i = 1,nmass
         if((qcflag.eq.1).or.(qmflag.eq.1)) then
            write(2,'(''MASS'',2x,i5,2x,'''',a4,2x,f8.5,2x,a4)')
     $           200+i,mass(i),setmas(mass(i)),mass(i)
         else 
            write(2,'(''MASS'',2x,i5,2x,''Y'',a4,2x,f12.8,2x,a4)')
     $           200+i,mass(i),setmas(mass(i)),mass(i)
         endif
      enddo
      write(2,*)
      write(2,'(''AUTO ANGL DIHE'')')
      write(2,*)
C
C     Find and store residue limits
C     CHECK: Do we ever need this... Right now this is not implemented!
C
      call reslim(natom,type,resnum,resinx)
C
C     Add hydrogens:
C     WARNING!!!!
C     The addhyd routine knows only a subset of atoms and their valencies!
C     It always assumes single bonds in the molecule, thus producing
C     the maximum number of hydrogens. It is up to user to delete the
C     excesive hydrogen atoms! Maybe HBUILD knows how to delete them??
C     It knows how to do for some of them!
C
C
c$$$C     TEST ADDATO():
c$$$      i=2320
c$$$      call addato('HYDR',i,natom,name,resn,resid,chain,type,
c$$$     $     x(i)+0.1,y(i),z(i),x,y,z)
c$$$C
C
      if(addhfl.eq.0) goto 89
      i=1
      cresid=resid(i)
      namgen=0
 88   continue
C
C     This puts unique atom names for added atoms per residue
      if (cresid.ne.resid(i)) then
         namgen=0
         cresid=resid(i)
      endif
C
      nadd = 0
      if((type(i).ne.'XXXX').and.(type(i)(1:1).ne.'H')) then
         numbon=0
         r1 = setrad(type(i))
         do j = i, natom
            if(type(j).eq.'XXXX')then
               r2 = setrad(name(j)(1:1)//'   ')
            else
               r2 = setrad(type(j))
            endif
            dist=dij(x(j),y(j),z(j),x(i),y(i),z(i))
C
C     TODO: double bonds, rings, triple bonds, ... ???
C           This is a good place to figure out double bonds???
C           One can get the distances for some of the pairs from
C           MOLDEN!!!
C
            if(dist.le.((r1+r2)/two)) then
               numbon=numbon+1
            endif
         enddo
         call addhyd(numbon,i,namgen,natom,nadd,type,
     $        x,y,z,name,resn,chain,resid)
      endif
C
      i=i+nadd+1
C
      if (natom.ge.maxat) then
         write(*,*)'Too many atoms: modify the program!'
         stop
      endif
C
      if (i.le.natom) goto 88
C
C     Update the chain info
      call seglim(ichain,nchain,chain,natom)
C
      write(*,*)'Number of atoms: ',natom
      write(*,*)'Number of chains/segments =',ichain
      write(*,*)'nchain(i)=',(nchain(i),i=1,ichain+1)
C
C     No addhyd needed:
 89   continue
C
C
C     Make residues: ATOM commands
C
C
CCC      open(unit=22,name='debug',status='unknown')
CCC      do i=1,natom
CCC         write(22,*)i,name(i),resid(i),resn(i)
CCC      enddo
      nato=0
      iresad=0
      resn(natom+1)='ZZZZ'
      ibres=1
      ieres=1
      newflag=.false.
      do i = 1, natom
         if (type(i).ne. 'XXXX') then
CCC            write(*,*)'nato,i,type,resn,resid,chain = ',
CCC     $           nato,i,type(i),resn(i),resid(i),chain(i)
C     In case of XYZ format we have charges at the end!
            ww=zero
            if (qxyz) ww=wmain(i)
            qnew=resinew(resn(i),resiad,iresad)
            if(qnew)newflag=.true.
            if (nato .eq. 0) then
               nato=nato+1
               if (newflag) then
CC                  write(2,'(''PATC FIRST NONE LAST NONE''/)')
                  write(2,'(''RESI'',2x,a4,2x,f3.1)')resn(i),0.0
                  if((qcflag.eq.1).or.(qmflag.eq.1)) then
                     write(2,'(''ATOM'',2x,a4,2x,'''',a4,2x,f7.2)')
     $                    name(i),type(i),ww
                  else 
                     write(2,'(''ATOM'',2x,a4,2x,''Y'',a4,2x,f7.2)')
     $                    name(i),type(i),ww
                  endif
                  ibres=i
                  ieres=i
C     If this is the only atom in the residue,
C     we can finish this residue
                  if (resid(i).ne.resid(i+1)) then
                     nato=0
                     write(2,*)
                  endif
               endif
            elseif (resid(i).eq.resid(i-1)) then
               if (newflag) then
                  ieres=ieres+1
                  if((qcflag.eq.1).or.(qmflag.eq.1)) then
                     write(2,'(''ATOM'',2x,a4,2x,'''',a4,2x,f7.2)')
     $                  name(i),type(i),ww
                  else 
                     write(2,'(''ATOM'',2x,a4,2x,''Y'',a4,2x,f7.2)')
     $                  name(i),type(i),ww
                  endif 
                  if(qxyz)then
                     xyzgrp=xyzgrp+1
                     if(xyzgrp.gt.80)then
                        write(2,'(''GROUP'')')
                        xyzgrp=0
                     endif
                  endif
                  if (resn(i).ne.resn(i+1)) then
                     nato=0
C     
C     Make BOND commands
C
                     if(qmol) then
                        do ix=1,nbonds
                           write(2,'(4a)')'BOND  ',
     $                          name(ibond(ix)),' ',name(jbond(ix))
                        enddo
                     else
                     do ix=ibres,ieres-1
                        ibflag=0
                        indx=0
CCC                        write(*,*)'ix,type=',ix,type(ix)
                        r1= setrad (type(ix))
                        do j=ix+1,ieres
                           dist = dij(x(ix),y(ix),z(ix),x(j),y(j),z(j))
CCC                           write(*,*)'j,type=',j,type(j)
                           r2= setrad (type(j))
                           requ =(r1+r2)/two
C
C     TODO:
C     Fix for hydrogens which are defined from experiment
C     For now we ignore them in the bonding ?!
C     0.51 comes from hand/(addhyd) editing of hydrogen atoms
C     for XYZ we treat hydrogen normally!
C
                           if ((.not.qxyz).and.
     $                          ((type(ix)(1:1).eq.'H').or.
     $                          (type(j)(1:1).eq.'H'))) requ = 0.51d0
CCC                           write(*,*)'ix,j,t(i),t(j),r1,r2,d=',ix,j,
CCC     $                       type(ix),type(j),r1,r2,dist
C     If both are hydrogens we don't want a bond there!
                           if ((type(ix)(1:1).eq.'H').and.
     $                          (type(j)(1:1).eq.'H')) requ = 0.1d0
                           if (dist .lt. requ) then
                              if(ibflag.eq.0)then
                                 ibflag = 1
                                 indx=0
                                 chafil(1:6) = 'BOND  '
                              end if
                              write(tmp,'(1x,a4,2x,a4,3x)')
     1                             name(ix), name(j)
                              chafil(7+indx*14:20+indx*14)=tmp(1:14)
                              indx=indx+1
                              if (indx.eq.5) then
                                 ibflag=0
                                 write(2,'(a)')chafil(1:76)
                              end if
                           end if
                        end do
                        indx = indx - 1
                        if(ibflag.eq.1) write(2,'(a)')
     $                       chafil(1:20+14*indx)
                     end do
                     endif
                     write(2,*)
                  endif
               else
                  if (resid(i).ne.resid(i+1)) nato=0
               endif
            endif
            if(resid(i).ne.resid(i+1))newflag=.false.
         endif
      enddo
C
C     End of BOND commands
C
      write(2,'(/''PATC FIRST NONE LAST NONE'')')
      write(2,*)
C
C     Make patches for bonds between residues of non 'XXXX' type,
C     if necessary
C
C     TODO: if HATATOM residue is of peptide type then one bond is
C           too much!!!!
C
C     TODO: Also when making this inter residue bonds, delete the hydrogens
C           if present !!! ???
C
C     CHECK:      This is good only for non 'H   ' types !!!!
C
      ipat=0
      do i = 1, ichain
         do j = nchain(i),nchain(i+1)-1
            if(type(j).eq.'XXXX')then
               r1 = setrad(name(j)(1:1)//'   ')
            else
               r1 = setrad(type(j))
            endif
            do k = j+1, nchain(i+1)
               if((resid(j).ne.resid(k)).and.
     $              ((type(j).ne.'H   ').and.(type(k).ne.'H   ')).and.
     $              ((type(j).ne.'XXXX').or.(type(k).ne.'XXXX'))) then
                  dist=dij(x(j),y(j),z(j),x(k),y(k),z(k))
                  if(type(k).eq.'XXXX')then
                     r2 = setrad(name(k)(1:1)//'   ')
                  else
                     r2 = setrad(type(k))
                  endif
CCC                  write(*,*)'j,k,t(j),t(k),r1,r2=',j,k,
CCC     $                 type(j),type(k),r1,r2
                  if(dist.le.((r1+r2)/two)) then
                     ipat=ipat+1
                     call letts(ipat,cc,lcc)
                     write(tmp,'(a,a,a)')
     $                    resn(j)(1:1),resn(k)(1:1),cc(1:lcc)
                     write(2,'(''pres '',a4,f5.2)')tmp,0.0
                     write(2,'(''bond 1'',a4,'' 2'',a4)')
     $                    name(j),name(k)
                     write(2,*)
                     write(pcomm(ipat),'(''patch '',a4,2x,
     $                    a4,1x,a4,2x,a4,1x,a4)')
     $                    tmp,chain(j),resid(j),chain(k),resid(k)
                  endif
               endif
            enddo
         enddo
      enddo
C
C     Finish RTF definitions
C
      write(2,*)
      write(2,'(/,''End'')')
      write(2,*)
      write(2,'(a,a)')'open read unit 1 card name ',parfile(1:lparf)
      write(2,'(a)')'read param card unit 1'
      write(2,*)
C
C     Make parameters: radii 
C                      and bonds between the HETATOM types and H
C
C-------------------------------------------------
C 
      write(2,*)
      write(2,'(a)')'read param card append'
      write(2,'(a)')'* Specify radia only for nonbond !'
      write(2,'(a)')'* bonds and angles are good for HBUILD only !'
      write(2,'(a)')'* It also avoids cutoff methods to be used'
      write(2,'(a)')'* while checking nuclear repulsion energy,etc...'
      write(2,'(a)')'*'
      write(2,*)
      write(2,'(a)') 'bonds ! Don''t believe these ones'
      if((qcflag.eq.1).or.(qmflag.eq.1)) then
         do i = 1,nmass
            do j = 1, nmass
               !write(2,'(''Y'',a4,''   YH  '',2f7.2)')mass(i),350.0,1.0
                  ! bonds for hydrogens
                  if (mass(i).eq.'H') then
                     if (mass(j).eq.'N') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.00
                     else if (mass(j).eq.'C') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.07
                     else if (mass(j).eq.'H') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,0.60
                     else if (mass(j).eq.'S') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.31
                     else if (mass(j).eq.'O') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,0.96
                     else if (mass(j).eq.'P') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.40
                     else if (mass(j).eq.'M') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.46
                     else
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.00
                     endif
                  ! bonds for nitrogen
                  else if (mass(i).eq.'N') then
                     if (mass(j).eq.'N') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.40
                     else if (mass(j).eq.'C') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.47
                     else if (mass(j).eq.'S') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.71
                     else if (mass(j).eq.'O') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.36
                     else if (mass(j).eq.'M') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.86
                     else if (mass(j).eq.'P') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.80
                     else if (mass(j).eq.'H') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.00
                     else
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.80
                     endif
                  ! bonds for carbon
                  else if (mass(i).eq.'C') then
                     if (mass(j).eq.'N') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.47
                     else if (mass(j).eq.'C') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.54
                     else if (mass(j).eq.'S') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.78
                     else if (mass(j).eq.'O') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.43
                     else if (mass(j).eq.'M') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.93
                     else if (mass(j).eq.'P') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.87
                     else if (mass(j).eq.'H') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.07
                     else
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.80
                     endif
                  ! bonds for sulfur
                  else if (mass(i).eq.'S') then
                     if (mass(j).eq.'N') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.71
                     else if (mass(j).eq.'C') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.78
                     else if (mass(j).eq.'S') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.02
                     else if (mass(j).eq.'O') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.67
                     else if (mass(j).eq.'M') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.17
                     else if (mass(j).eq.'P') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.11
                     else
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.10
                     endif
                  ! bonds for oxygen
                  else if (mass(i).eq.'O') then
                     if (mass(j).eq.'N') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.36
                     else if (mass(j).eq.'C') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.43
                     else if (mass(j).eq.'S') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.67
                     else if (mass(j).eq.'O') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.32
                     else if (mass(j).eq.'M') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.82
                     else if (mass(j).eq.'P') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.76
                     else if (mass(j).eq.'H') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,0.96
                     else
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.60
                     endif
                  ! bonds for Phosphorus
                  else if (mass(i).eq.'P') then
                     if (mass(j).eq.'N') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.80
                     else if (mass(j).eq.'C') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.87
                     else if (mass(j).eq.'S') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.11
                     else if (mass(j).eq.'O') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.76
                     else if (mass(j).eq.'M') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.26
                     else if (mass(j).eq.'P') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.20
                     else if (mass(j).eq.'H') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.40
                     else
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.10
                     endif
                  else
                  ! bonds for Nitrogen
                     if (mass(j).eq.'N') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.70
                     else if (mass(j).eq.'C') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.80
                     else if (mass(j).eq.'H') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.00
                     else if (mass(j).eq.'S') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.10
                     else if (mass(j).eq.'O') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.70
                     else if (mass(j).eq.'M') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.20
                     else if (mass(j).eq.'P') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.10
                     else if (mass(j).eq.'H') then
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.00
                     else
                        write(2,'('''',a4,'''',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.10
                     endif
                  endif
            enddo
         enddo
      else
         do i = 1,nmass
            do j = 1, nmass
               !write(2,'(''Y'',a4,''   YH  '',2f7.2)')mass(i),350.0,1.0
                  ! bonds for hydrogens
                  if (mass(i).eq.'H') then
                     if (mass(j).eq.'N') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.00
                     else if (mass(j).eq.'C') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.07
                     else if (mass(j).eq.'H') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,0.60
                     else if (mass(j).eq.'S') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.31
                     else if (mass(j).eq.'O') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,0.96
                     else if (mass(j).eq.'P') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.40
                     else if (mass(j).eq.'M') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.46
                     else
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.00
                     endif
                  ! bonds for nitrogen
                  else if (mass(i).eq.'N') then
                     if (mass(j).eq.'N') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.40
                     else if (mass(j).eq.'C') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.47
                     else if (mass(j).eq.'S') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.71
                     else if (mass(j).eq.'O') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.36
                     else if (mass(j).eq.'M') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.86
                     else if (mass(j).eq.'P') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.80
                     else if (mass(j).eq.'H') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.00
                     else 
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.80
                     endif 
                  ! bonds for carbon
                  else if (mass(i).eq.'C') then
                     if (mass(j).eq.'N') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.47
                     else if (mass(j).eq.'C') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.54
                     else if (mass(j).eq.'S') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.78
                     else if (mass(j).eq.'O') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.43
                     else if (mass(j).eq.'M') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.93
                     else if (mass(j).eq.'P') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.87
                     else if (mass(j).eq.'H') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.07
                     else 
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.80
                     endif
                  ! bonds for sulfur 
                  else if (mass(i).eq.'S') then
                     if (mass(j).eq.'N') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.71
                     else if (mass(j).eq.'C') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.78
                     else if (mass(j).eq.'S') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.02
                     else if (mass(j).eq.'O') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.67
                     else if (mass(j).eq.'M') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.17
                     else if (mass(j).eq.'P') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.11
                     else
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.10
                     endif
                  ! bonds for oxygen
                  else if (mass(i).eq.'O') then
                     if (mass(j).eq.'N') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.36
                     else if (mass(j).eq.'C') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.43
                     else if (mass(j).eq.'S') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.67
                     else if (mass(j).eq.'O') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.32
                     else if (mass(j).eq.'M') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.82
                     else if (mass(j).eq.'P') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.76
                     else if (mass(j).eq.'H') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,0.96
                     else
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.60
                     endif
                  ! bonds for Phosphorus
                  else if (mass(i).eq.'P') then
                     if (mass(j).eq.'N') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.80
                     else if (mass(j).eq.'C') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.87
                     else if (mass(j).eq.'S') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.11
                     else if (mass(j).eq.'O') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.76
                     else if (mass(j).eq.'M') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.26
                     else if (mass(j).eq.'P') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.20
                     else if (mass(j).eq.'H') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.40
                     else
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.10
                     endif
                  else 
                  ! bonds for Nitrogen
                     if (mass(j).eq.'N') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.70
                     else if (mass(j).eq.'C') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.80
                     else if (mass(j).eq.'H') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.00
                     else if (mass(j).eq.'S') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.10
                     else if (mass(j).eq.'O') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.70
                     else if (mass(j).eq.'M') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.20
                     else if (mass(j).eq.'P') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.10
                     else if (mass(j).eq.'H') then
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,1.00
                     else
                        write(2,'(''Y'',a4,''Y'',a4,2f7.2)')
     &                  mass(i),mass(j),350.0,2.10
                     endif
                  endif
            enddo 
         enddo
      endif
C     
C     Add also angle parameters:
C
      write(2,*)
      write(2,'(a)') 'angles ! Don''t believe these ones'
      if((qcflag.eq.1).or.(qmflag.eq.1)) then
         do i=1,nmass
            do j=1,nmass
               do k=1,nmass
                  write(2,'('''',a4,3x,'''',a4,3X,'''',a4,3x,2f7.2)')
     $                  mass(i),mass(j),mass(k), 50.0, 110.0
               enddo 
            enddo 
         enddo 
!        do i = 1, nmass
!           do j = 1, nmass
!              write(2,'('''',a4,3x,'''',a4,''   H  '',2f7.2)')
!    $              mass(i),mass(j), 50.0, 110.0
!           enddo  
!        enddo
      else 
         do i=1,nmass
            do j=1,nmass
               do k=1,nmass
                  write(2,'(''Y'',a4,3x,''Y'',a4,3X,''Y'',a4,3x,2f7.2)')
     $                  mass(i),mass(j),mass(k), 50.0, 110.0
               enddo 
            enddo 
         enddo 
!        do i = 1, nmass
!           do j = 1, nmass
!              write(2,'(''Y'',a4,3x,''Y'',a4,''   YH  '',2f7.2)')
!    $              mass(i),mass(j), 50.0, 110.0
!           enddo  
!        enddo
      endif 
C     
C     Add some dihedral parameters:
C
      write(2,*)
      write(2,'(a)') 'dihedral ! Don''t believe these ones'
      if((qcflag.eq.1).or.(qmflag.eq.1)) then
         do i=1,nmass
            do j=1,nmass
            write(2,'(''X'',7x,'''',a4,3X,'''',a4,3x,''X'',3x,2f7.2)')
     $                  mass(i),mass(j), 5.0 1,  1.0
            enddo
         enddo
      else 
         do i=1,nmass
            do j=1,nmass
            write(2,'(''X'',7x,''Y'',a4,3X,''Y'',a4,3x,''X'',3x,2f7.2)')
     $                  mass(i),mass(j), 5.0 1,  1.0
            enddo
         enddo

      endif
C     
      write(2,*)
      write(2,'(a)')'NBOND ! CUTNB 999.0 CTONNB 990.0 CTOFNB 995.0'
C
      if((qcflag.eq.1).or.(qmflag.eq.1)) then
         do i = 1,nmass
            write(2,'('''',a4,3f8.5)')mass(i),0.0,0.0, setrad(mass(i))
         enddo
      else 
         do i = 1,nmass
            write(2,'(''Y'',a4,3f8.5)')mass(i),0.0,0.0, setrad(mass(i))
         enddo
      endif 
      write(2,*)
      write(2,'(a)')'End'
C----------------------------------------------------
C
C     Write a .CRD file, and parse for sequence(s)
C
      if(.not.qxyz) then
         write(4,'(''* '',a,
     $        '': automagically generated from PDB,XYZ,or other format''
     $        )')pdb
         write(4,'(''*'')')
      endif
      write(4,'(i5)')natom
      ires=1
      cresn=resn(1)
      cresid=resid(1)
      seq(ires)=resn(1)
      cseg=chain(1)
      iseg=1
      lseq(iseg)=ires
      nseg(1)=cseg
      do i = 1, natom
         if ((resn(i).ne.cresn).or.(resid(i).ne.cresid)) then
            cresn=resn(i)
            cresid=resid(i)
            ires=ires+1
            seq(ires)=resn(i)
         endif
         if (chain(i).ne.cseg) then
            cseg=chain(i)
            iseg=iseg+1
            lseq(iseg)=ires
            nseg(iseg)=cseg
         endif
         write(4,'(i5,i5,1x,a4,1x,a4,3f10.5,1x,a4,1x,a4,f10.5)')
     $        i, ires, resn(i), name(i), x(i), y(i), z(i), chain(i),
     $        resid(i), wmain(i)
      end do
C 
C-------------------------------------------------
C
C     Make PSF (via read sequ & generate commands)
C
      lseq(iseg+1)=ires+1
CCC      write(*,*)'Number of segments = ',iseg
      do i=1,iseg
CCC         write(*,*)'i,lseq,ires=',i,lseq(i),ires
         write(2,*)
         if (nseg(i)(1:2).eq.'CW') then
            write(2,'(''!read sequ tip3 '',i5)') lseq(i+1)-lseq(i)
            write(2,*)
            call letts(i,cc,lcc)
            write(tmp,'(a,''-'',a,''.PDB'')')filnam(1:basenm),cc(1:lcc)
            write(2,'(''open read unit 1 card name "'',a,''"'')')
     $           tmp(1:basenm+5+lcc)
            write(2,'(''read sequ pdb unit 1'')')
            write(2,*)
            write(2,'(''generate '',a4,'' noangle nodihe'')') nseg(i)
         else
            write(2,'(''!read sequ card'')')
            write(2,'(''!*'')')
            write(2,'(''!'',i4)') lseq(i+1)-lseq(i)
            write(2,'(''!'',20a4)')(seq(k), k=lseq(i),lseq(i+1)-1)
            write(2,*)
            call letts(i,cc,lcc)
            write(tmp,'(a,''-'',a,''.PDB'')')filnam(1:basenm),cc(1:lcc)
            write(2,'(''open read unit 1 card name "'',a,''"'')')
     $           tmp(1:basenm+5+lcc)
            write(2,'(''read sequ pdb unit 1'')')
            write(2,*)
            write(2,'(''generate '',a4,'' setup'')') nseg(i)
         endif
         write(2,*)
      enddo


ccccccccccccccccccccccccccccccccccccc
c      Add S-S bonds
ccccccccccccccccccccccccccccccccccccc
      
      do i=1,natom-1
         if ((type(i).eq.'XXXX').and.(name(i).eq.'SG  ')) then
            r1= setrad ('S   ')
            do j=i+1,natom
c               r2= setrad (type(j))
               if ((type(j).eq.'XXXX').and.(name(j).eq.'SG  ')) then
                  r2= setrad ('S   ')

                  dist = dij(x(i),y(i),z(i),x(j),y(j),z(j))
c     write(*,*) dist
                  requ =(r1+r2)/two
c                  write(*,*) 'ti,tj,ri,rj,r=',type(i),type(j),r1,r2,requ
                  if (dist .lt. 2.5) then
                     write(2,'(''patch disu'',1x,a4,1x,a4,1x,a4,1x,a4)')
     $                    chain(i),resid(i),chain(j),resid(j)
                     
                  endif
               endif
            enddo
         endif
      enddo
      write(2,*)
C
C     Add patches created above, if any
C
      if (ipat.gt.0) then
         do i = 1, ipat
            write(2,'(a)')pcomm(i)
         enddo
      endif
C
C     Read the created coordinate (.CRD) file
C
      jlen=islen(xname)
      write(2,*)
      write(2,'(''open read unit 1 card name "'',a,''"'')')xname(1:jlen)
      write(2,'(a)')'read coor card unit 1'
      write(2,*)
C
C     TODO:
C     top_all22_prot file has a BUG in C-terminus oxygen IC table
C     We can fix it with the IC EDIT command
C
C     Build the atoms, which are not available from experiment
C
      write(2,'(a)')'ic param'
      write(2,'(a)')'ic build'
      write(2,*)
C
C     Clear the coordinates of all the hydrogens, and use HBUILD
C     to make new hydrogens.
C     TODO:
C           This procedure should be protected by command line
C           parameter so it does not build the hydrogen in the
C           cases where they were assigned by some other method.
C
      write(2,'(a)')'coor init sele type h* end'
      write(2,'(a)')'hbuild'
      write(2,*)
C
C     Clean the molecule: This also shows how many atoms didn't
C     make it into molecule !?
C
      write(2,'(a)')'delete atom sele .not. init show end'
      write(2,*)
C
C     Store the resulting coordnates, which includes all the hydrogens
C
C     write(2,'(''open write unit 1 card name '',a)') xname
C     write(2,'(a)')'write coor card unit 1'
C     write(2,*)

C
C     Tim Miller (6/30/2008): Rename the seg if desired
C
      if(rnmeseg.eqv..true.) then
         write(2,'(''rename segid '',a,'' select all end'')')segname
      endif

C
C     Store the resulting coordinates, which includes all the hydrogens
C
      write(2,'(''open write unit 1 card name '',a,''-f.crd'')')
     $     xname(1:basenm)
      write(2,'(a)')'write coor card unit 1'
      write(2,*)

C
C     Store the resulting PSF, which includes all the hydrogens
C
      write(2,'(''open write unit 1 card name '',a,''-f.psf'')')
     $     xname(1:basenm)
      write(2,'(a)')'write psf card unit 1'
      write(2,*)
C
C     Store the resulting PDB.
C
      write(2,'(''open write unit 1 card name '',a,''-f.pdb'')')
     $     xname(1:basenm)
      write(2,'(a)')'write coor pdb unit 1'
      write(2,*)
C
C     Store the resulting coordinates in the trajectory format for VMD
C
      write(2,
     $     '(''! See these files with: vmd -psf hb'',a,
     $     ''.psf -dcd hb'',a,''.dcd'')')xname(1:basenm),xname(1:basenm)
      write(2,'(''open write unit 1 file name hb'',a,''.dcd'')')
     $     xname(1:basenm)
      write(2,'(a)')'traj iwrite 1 nwrite 1 nfile 1'
      write(2,'(a)')'traj write'
      write(2,*)
C  
CCC      write(*,*)'i,lseq=',iseg+1,lseq(iseg+1)
      if (qmflag.eq.1) then
         write(2,'(3a)')'envi input   "',filnam(1:basenm),'.str"'
         write(2,'(3a)')'envi output  "',filnam(1:basenm),'.gms"'
         write(2,'(3a)')'envi punch   "',filnam(1:basenm),'.dat"'
         write(2,'(a)')'envi dictnry "test.f10"'
         write(2,'(a)')'envi work15  "test.f15"'
         write(2,'(a)')'envi dasort  "test.f20"'
         write(2,*)
         write(2,'(a)')'gamess remove exgr noguess select all end'
         write(2,*)
         write(2,'(a)')'energy'
         write(2,*)
         open(unit=20,file=filnam(1:basenm)//'.str',status='unknown')
         write(20,'(a,i3)')' $CONTRL COORD=CART NOSYM=1 ICHARG=',icharg
         write(20,'(a)')'          SCFTYP=RHF MAXIT=50'
         write(20,'(a)')' !         runtyp=prop'
         write(20,'(a)')' !         nprint=-5'
         write(20,'(a)')' !         RUNTYP=GRADIENT'
         write(20,'(a)')' !         RUNTYP=OPTIMIZE'
         write(20,'(a)')' !         EXETYP=CHECK'
         write(20,'(a)')' !         MOLPLT=.TRUE.'
         write(20,'(a)')' $END'
         write(20,'(a)')' $SYSTEM MEMORY=2000000 TIMLIM=1000000 $END'
         write(20,'(a)')' $BASIS'
         write(20,'(a)')' !          GBASIS=AM1'
         write(20,'(a)')' !          GBASIS=STO NGAUSS=3'
         write(20,'(a)')'           GBASIS=N31 NGAUSS=6 '
         write(20,'(a)')'           NDFUNC=1'
         write(20,'(a)')' !          NDFUNC=1 NPFUNC=1'
         write(20,'(a)')' !          DIFFSP=.TRUE. DIFFS=.TRUE.'
         write(20,'(a)')' $END'
         write(20,'(a)')' $SCF  DIRSCF=.True. soscf=.true. $END'
         write(20,'(a)')' ! $STATPT  NSTEP=100 OPTTOL=0.0004 $END'
         write(20,'(a)')' $DATA'
         write(20,*)
         write(20,*)
         write(20,*)
         write(20,'(a)')' $END'
         write(20,'(a)')' $ELPOT IEPOT=1 where=pdc output=punch $END'
         write(20,'(a)')' $PDC CONSTR=NONE $END'
         close(20)
      endif
      
      if (qcflag.eq.1) then
         write(2,'(3a)')'envi qchemcnt "qchem.inp"'
         write(2,'(3a)')'envi qchemexe "qchem"'
         write(2,'(3a)')'envi qchemout "q1.out"'
         write(2,'(3a)')'envi qcheminp "q1.inp"'
         write(2,'(3a)')'! Need this for Q-Chem Replica/Path'
         write(2,'(3a)')'envi qchempwd "current directory"'
         write(2,'(3a)')'envi PBS_NODEFILE "qchosts"'
         write(2,'(3a)')' '
         write(2,*)' ' 
         write(2,'(a)')'qchem remove select all end'
         write(2,*)' ' 
         write(2,'(a)')'energy'
         write(2,*)' '
         write(2,'(a)')'stop'
         open(unit=21,file='qchem.inp',status='unknown')
         write(21,'(a)')'$comment'
         write(21,'(a)')'qchem control file needed.'
         write(21,'(a)')'$end'
         write(21,'(a)')' '
         write(21,'(a)')'$rem'
         write(21,'(a)')'exchange              hf'
         write(21,'(a)')'basis                 6-31G*'
         write(21,'(a)')'qm_mm                 true'
         write(21,'(a)')'jobtype               force'
         write(21,'(a)')'symmetry              off'
         write(21,'(a)')'sym_ignore            true'
         write(21,'(a)')'print_input           false'
         write(21,'(a)')'qmmm_print            true'
         write(21,'(a)')'scf_convergence       6'
         write(21,'(a)')'scf_algorithm         diis'
         write(21,'(a)')'MAX_SCF_CYCLES        1000'
         write(21,'(a)')'MAXKL                 150'
         write(21,'(a)')'$end'
         write(21,'(a)')' '
         write(21,'(a)')'$molecule'
         write(21,'(a,i3)')'0 ',icharg
         write(21,'(a)')'$end'
         write(21,'(a)')' '
         write(21,'(a)')' '
         write(21,'(a)')' '
         close(21)
      endif
      
      stop
C
      end
C
      subroutine clrstr(n,a)
      implicit none
      character*4 a(*)
      integer i,n
      do i=1,n
         a(i)='    '
      enddo
      return
      end

      real*8 function dij(x1,y1,z1,x2,y2,z2)
      implicit none
      real*8 x1,y1,z1,x2,y2,z2
      real*8 s1,s2,s3
      s1=x1-x2
      s2=y1-y2
      s3=z1-z2
c      write(*,*) s1,s2,s3
      dij = sqrt(s1*s1+s2*s2+s3*s3)
      return
      end

      integer function islen(a)
      implicit none
      character*(*) a
      integer i,k
      i=0
 1    i=i+1
      k=ichar(a(i:i))
C     write(*,*)'k,i=',k,i
      if ( (k.le.32) .or. (k.gt.127) ) goto 2
      goto 1
 2    islen = i - 1
      return
      end

      subroutine uc(src,dst)
      implicit none
      character*(*) src,dst
      integer len,i,islen
      len=islen(src)
      do i = 1, len
         dst(i:i)=src(i:i)
         if (src(i:i).ge.'a') dst(i:i)=char(ichar(src(i:i))-32)
      end do
      if (len.lt.4) dst(len+1:4)=' '
      return
      end
C
      subroutine strcln(id)
      implicit none
      character*4 id, tmp
      integer i,j
C
      tmp(1:4)=id(1:4)
      id(1:4)='    '
      j=1
      do i =1,4
         if(tmp(i:i).ne.' ') then
            id(j:j)=tmp(i:i)
            j=j+1
         endif
      enddo
      return
      end
C
      logical function isdigit(c)
      implicit none
      character*1 c
      isdigit=.false.
      if((ichar(c).ge.ichar('0')).and.(ichar(c).le.ichar('9')))
     $     isdigit=.true.
      return
      end
C
      subroutine namcln(id)
      implicit none
      character*4 id
      character*1 t
      integer i,j,linlen
      logical isdigit
C
C     In PDB sometimes we have like 1H6 names
C     Change this to H61
C
      if(isdigit(id(1:1)))then
         i=linlen(id)
         t=id(1:1)
         do j=1,i-1
            id(j:j)=id(j+1:j+1)
         enddo
         id(i:i)=t
      endif
C
      return
      end
C
      subroutine strlj(id)
      implicit none
      character*4 id, tmp
      integer i,j
C
      tmp(1:4)=id(1:4)
      id(1:4)='    '
      j=4
      do i =4,1,-1
         if(tmp(i:i).ne.' ') then
            id(j:j)=tmp(i:i)
            j=j-1
         endif
      enddo
      return
      end
C
      subroutine terfix(natom,name,resid,type,chain)
      implicit none
      integer natom,i,j,k
      character*4 name(*),resid(*),type(*),chain(*)
C
      do i = 1, natom
         if(name(i).eq.'OXT ') then
            j=i-1
            k=i+1
            do while (resid(j).eq.resid(i))
               if (name(j).eq.'O   ') then
                  name(j)='OT2 '
               endif
               j=j-1
            enddo
            do while (resid(k).eq.resid(i))
               if (name(k).eq.'O   ') then
                  name(k)='OT2 '
               endif
               k=k+1
            enddo
            name(i)='OT1 '
         endif
      enddo
C
      do i=1,natom
         if((chain(i).ne.chain(i+1)).and.(type(i).eq.'XXXX')) then
            k=i-1
            do while (resid(k).eq.resid(i))
               if (name(k).eq.'O   ') then
                  name(k)='OT2 '
               endif
               k=k-1
            enddo
         endif
      enddo
      return
      end


      subroutine settyp(mode,name,type)
      implicit none
      character*4 name,type
      integer mode
      if (mode.eq.0) then
         type='XXXX'
         return
      else if (mode.eq.1) then
         type(1:4)='    '
         type(1:2)=name(1:2)
         return
      else if (mode.eq.2) then
         type=name(1:2)
         return
      endif
      return
      end

      subroutine atonam(atoms)
      implicit none
      integer maxper,i
      parameter (maxper=103)
      character*4 atoms(*),atom(maxper)
C
      data atom/
     $   'H   ','HE  ','LI  ','BE  ','B   ','C   ','N   ','O   ','F   ',
     $   'NE  ','NA  ','MG  ','AL  ','SI  ','P   ','S   ','CL  ','AR  ',
     $   'K   ','CA  ','SC  ','TI  ','V   ','CR  ','MN  ','FE  ','CO  ',
     $   'NI  ','CU  ','ZN  ','GA  ','GE  ','AS  ','SE  ','BR  ','KR  ',
     $   'RB  ','SR  ','Y   ','ZR  ','NB  ','MO  ','TC  ','RU  ','RH  ',
     $   'PD  ','AG  ','CD  ','IN  ','SN  ','SB  ','TE  ','I   ','XE  ',
     $   'CS  ','BA  ','LA  ','CE  ','PR  ','ND  ','PM  ','SM  ','EU  ',
     $   'GD  ','TB  ','DY  ','HO  ','ER  ','TM  ','YB  ','LU  ','HF  ',
     $   'TA  ','W   ','RE  ','OS  ','IR  ','PT  ','AU  ','HG  ','TL  ',
     $   'PB  ','BI  ','PO  ','AT  ','RN  ','FR  ','RA  ','AC  ','TH  ', 
     $   'PA  ','U   ','NP  ','PU  ','AM  ','CM  ','BK  ','CF  ','ES  ',
     $   'FM  ','MD  ','NO  ','LR  '/        
C
      do i = 1, maxper
         atoms(i)=atom(i)
      enddo
      return
      end

      real*8 function setmas(type)
      implicit none
      integer i,j,maxper
      parameter (maxper=103)
      character*4 type,atoms(maxper)
      real*8 masses(maxper),zero
C
      data masses/
     $     1.00794   ,4.002602  ,6.941    ,9.01218  ,10.811   ,
     $     12.011    ,14.00674  ,15.9994  ,18.998403,20.1797  ,
     $     22.989768 ,24.305    ,26.981539,28.0855  ,30.973762,
     $     32.066    ,35.4527   ,39.948   ,39.0983  ,40.078   ,
     $     44.95591  ,47.88     ,50.9415  ,51.9961  ,54.93805 ,
     $     55.847    ,58.9332   ,58.6934  ,63.546   ,65.39    ,
     $     69.723    ,72.61     ,74.92159 ,78.96    ,79.904   ,  
     $     83.8      ,85.4678   ,87.62    ,88.90585 ,91.224   ,
     $     92.90638  ,95.94     ,97.9072  ,101.07   ,102.9055 ,
     $     106.42    ,107.8682  ,112.411  ,114.818  ,118.71   ,
     $     121.760   ,127.6     ,126.90447,131.29   ,132.90543,
     $     137.327   ,138.9055  ,140.115  ,140.90765,144.24   ,
     $     144.9127  ,150.36    ,151.965  ,157.25   ,158.92534,
     $     162.50    ,164.93032 ,167.26   ,168.93421,173.04   ,
     $     174.967   ,178.49    ,180.9479 ,183.84   ,186.207  ,
     $     190.23    ,192.22    ,195.08   ,196.96654,200.59   ,
     $     204.3833  ,207.2     ,208.98037,208.9824 ,209.9871 ,
     $     222.0176  ,223.0197  ,226.0254 ,227.0278 ,232.0381 ,
     $     231.03588 ,238.0289  ,237.048  ,244.0642 ,243.0614 ,
     $     247.0703  ,247.0703  ,251.0796 ,252.083  ,257.0951 ,
     $     258.1     ,259.1009  ,262.11  /   
      data zero /0.0d0/
C
      call atonam(atoms)
C
      do j=1,9
         if ((ichar(type(2:2)) -(ichar("0"))).eq. j) then
            type=type(1:1)
         else
            type=type(1:2)
         endif 
      enddo

      setmas=zero
      do i=1,maxper
        if(type(1:1) .eq. atoms(i)) setmas = masses(i)
      enddo
      return
      end


      real*8 function setrad(type)
      implicit none
      integer maxper
      parameter (maxper=103)
      character*4 type,atoms(maxper)
      integer i
C
C     TODO: Find radii for upto maxper
      real*8 waals(54)
C
C     nanotube correction setrad('C')=1.75 (was 1.7)
C
      data waals /
     $     1.20d0,                                          1.40d0,
     $     1.82d0,1.78d0,1.74d0,1.90d0,1.55d0,1.52d0,1.47d0,1.54d0,
     $     2.27d0,2.22d0,2.16d0,2.10d0,1.80d0,2.00d0,1.75d0,1.88d0,
     $     2.75d0,2.57d0,
     $     2.56d0,2.54d0,2.52d0,2.50d0,2.48d0,2.46d0,2.44d0,2.42d0,
     $     2.41d0,2.40d0,
     $     2.40d0,2.10d0,1.85d0,1.90d0,1.85d0,2.02d0,
     $     3.10d0,2.80d0,
     $     2.77d0,2.74d0,2.71d0,2.68d0,2.65d0,2.62d0,2.59d0,2.56d0,
     $     2.53d0,2.51d0,
     $     2.50d0,2.20d0,2.10d0,2.06d0,1.98d0,2.16d0 /
C
      if (type.eq.'XXXX') then
         write(*,*)'This type is not supported!'
         write(*,*)'Type = ',type
         stop 'Stopped by wrong type error'
      endif
C
      call atonam(atoms)
C
C     Default radius is 2.0 A
C
      setrad=2.5d0
C
      do i=1,54
         if (type .eq. atoms(i)) setrad = waals(i)
      enddo
      return
      end
C
      subroutine addhyd(nb,i,namgen,natom,nadd,
     $     type,x,y,z,name,resn,chain,resid)
      implicit none
C
C     Routine which adds hydrogens to hetatoms. It assumes
C     there are only single bonds in the residue, but maybe we can
C     later delete the excessive hydrogen by some automatic procedure???
C
C     nb    - number of bonds
C     i     - current index in the arrays
C     natom - number of atoms: this is incremented by the number
C             of hydrogens found in this routine
C     the rest of parameters are from 1 .. natom
C
      integer nb,i,namgen,natom,nadd
      real*8 x(*),y(*),z(*)
      character*4 type(*),name(*),resn(*),chain(*),resid(*)
C
      integer j
      character*4 nat
C
      nadd=0
      if(type(i).eq.'XXXX') return
      if(nb.eq.0) return
C
C     Find the (maximum) number of hydrogens for HETATOM i
C
      nadd=10
C
      if(type(i)(1:1).eq.'C') nadd=4-nb
C
C     valence for our Nitrogen is 3 !?
      if(type(i)(1:1).eq.'N') nadd=3-nb
C
      if(type(i)(1:1).eq.'O') nadd=2-nb
C
C     Valence for Sulfur is 2 !?
      if(type(i)(1:1).eq.'S') nadd=2-nb
C
C     Valence for Phosphorus is 4 !?
      if(type(i)(1:1).eq.'P') nadd=4-nb
C
C     No hydrogen on sodium:
      if(type(i)(1:2).eq.'NA') nadd=0
C
C     This is for the existing Hydrogens
      if(type(i)(1:1).eq.'H') nadd=1-nb
C
C     The following line protects for cases were more
C     bonds were found then the atom valence.
C     Set nadd to zero in this case:
      if(nadd.lt.0)nadd=0
C
      if(nadd.eq.0) return
      if(nadd.eq.10) then
         nadd=0
         write(*,*)'Unknown type in addhyd for  HETATOM!'
         write(*,*)'Modify the source to correct this'
         write(*,*)'For now we ignore the routine and continue'
         write(*,*)'nb,i,name,type=',nb,i,name(i),type(i)
         return
      endif
C
C     add nadd atoms:
C
      do j = 1, nadd
         nat(1:1)='H'
         nat(2:3)=name(i)(1:2)
         if(namgen.lt.26) then
            nat(4:4)=char(ichar('A')+namgen)
         else
            nat(3:3)=char(ichar('A')+namgen/25)
            nat(4:4)=char(ichar('A')+mod(namgen,25))
         endif
         call strcln(nat)
         call addato(nat,i,natom,name,resn,resid,chain,type,
     $        x(i)+j*(0.5/nadd),y(i),z(i),x,y,z)
         namgen=namgen+1
      enddo
C
      return
      end
C
      subroutine letts(n,c,l)
      implicit none
      integer n,l,m,k
      character*(*) c
C
C     For now this works only for 3 characters, but the call
C     is ready for general number. Only the implementation needs
C     to be extended if needed, of course
C
      m=3
      if (n.gt.26**m) then
         write(*,*)'Too big number for this name: ',n
         stop
      endif
      l=1
      if(n.le.26) then
         c(1:1)=char(ichar('A')+n-1)
      else if (n.le.26**2+26) then
         c(1:1)=char(ichar('A')+(n-1)/26-1)
         c(2:2)=char(ichar('A')+mod((n-1),26))
         l=2
      else
         k=(n-1)/(26**2+26)-1
         c(1:1)=char(ichar('A')+k)
         k=mod(n-1,26**2+26)
         k=(k-1)/27
         c(2:2)=char(ichar('A')+k)
         k=mod(n-1,26)
         c(3:3)=char(ichar('A')+k)
         l=3
      endif
C
      return
      end
C
      subroutine reslim(natom,type,resnum,resinx)
      implicit none
      integer maxat
      parameter (maxat=20000)
C
C     CHECK: Do we really need this routine?
C            It is empty for now!!!
C
C     This routine finds the limits of the residues
C     for the whole system
C
      integer natom,resnum,resinx(*)
c     character*4 type
      character*4 type(maxat)

C
      return
      end
C
      integer function linlen(a)
C
C     This routine finds out the length of C*(*) a
C     without the blanks at the end of it. All the
C     blanks in the middle of it are not ignored!
C
      implicit none
      character*(*) a
      integer i,k
      i=len(a)+1
 1    i=i-1
      if (i.le.0) goto 2
      k=ichar(a(i:i))
C     write(*,*)'k,i=',k,i
      if (k.eq.32)  goto 1
 2    continue
      linlen = i
      return
      end
C
      subroutine addato(atn,icurr,natom,name,resn,resid,chain,type,
     $     cx,cy,cz,x,y,z)
      implicit none
C
C     This routine adds an atom after the icurr-th atom
C     It has a name atn
C     Inherits resn,resid,chain
C     It has to update the following parameters:
C           natom,name,resn,resid,chain,x,y,z,type
C
      integer icurr,natom
      real*8 cx,cy,cz,x(*),y(*),z(*)
      character*4 atn,name(*),resn(*),resid(*),chain(*),type(*)
C
      integer i
C
C     Make space for one atom
C
      natom=natom+1
C
C     chain() is special for ending info!
C
      chain(natom+1)=chain(natom)
C
      do i = natom,icurr+2,-1
         name(i)=name(i-1)
         resn(i)=resn(i-1)
         resid(i)=resid(i-1)
         chain(i)=chain(i-1)
         type(i)=type(i-1)
         x(i)=x(i-1)
         y(i)=y(i-1)
         z(i)=z(i-1)
      enddo
C
C     Add atn into its own position
C
      name(icurr+1)=atn
      resn(icurr+1)=resn(icurr)
      resid(icurr+1)=resid(icurr)
      chain(icurr+1)=chain(icurr)
      call settyp(1,atn,type(icurr+1))
      x(icurr+1)=cx
      y(icurr+1)=cy
      z(icurr+1)=cz
C
      return
      end
C
      subroutine seglim(ichain,nchain,chain,natom)
      implicit none
C
      integer ichain,natom,nchain(*)
      character*4 chain(*)
C
      integer i
C
      ichain=1
      nchain(ichain)=1
      do i = 2, natom
         if (chain(i).ne.chain(i-1)) then
            ichain=ichain+1
            nchain(ichain)=i
         endif
      enddo
      nchain(ichain+1)=natom
C
      return
      end
C
      subroutine chkchain(k,ichain,nchain,chain,newch)
      implicit none
C
      integer k, ichain, nchain(*)
      character*(*) chain(*),newch(*)
C
      integer i,linlen,l,lcc
      character*4 tmp
      character*2 cc
C
      l=linlen(chain(k))
      newch(k)=chain(k)
      if(ichain.lt.2) return
      tmp(1:l)=chain(k)(1:l)
CCC      write(*,*)'chkchain>k,ichain,tmp=',k,ichain,tmp(1:l)
      do i=1,ichain-1
         if(chain(nchain(i))(1:l).eq.tmp(1:l))then
            call letts(ichain,cc,lcc)
            write(newch(k),'(a,a)')tmp(1:l),cc(1:lcc)
CC            write(*,*)'chkchain>k,ichain,tmp,chain=',
CC     $           k,ichain,tmp(1:l),chain(k)(1:l+lcc)
            return
         endif
      enddo
C
      return
      end
C
      logical function resinew(resn,resiad,k)
      implicit none
C
      character*(*) resn
      character*(*) resiad(*)
      integer k,i
C
      if(k.gt.0) then
         do i = 1, k
            if (resn(1:4).eq.resiad(i)(1:4)) then
               resinew=.false.
               return
            endif
         enddo
      endif
      k=k+1
      resiad(k)(1:4)=resn(1:4)
      resinew=.true.
C
      return
      end
C
      logical function topres(resn)
      implicit none
C
      integer maxres
      parameter(maxres=30)
      character*4 resn,top(maxres)
      integer i
      logical t
C
C     TODO: add also stuff for DNA etc...
C
      data top/'ALA ', 'ARG ', 'ASN ', 'ASP ', 'CYS ', 'GLN ', 'GLU ',
     $         'GLY ', 'HSD ', 'HSE ', 'HSP ', 'ILE ', 'LEU ', 'LYS ',
     $         'MET ', 'PHE ', 'PRO ', 'SER ', 'THR ', 'TRP ', 'TYR ',
     $         'VAL ', 'TIP3', 'HOH ', 'CAL ', 'HIS ', 'A   ', 'C   ',
     $         'G   ', 'T   '/
C
      t=.false.
      do i=1,maxres
         if (resn.eq.top(i))t=.true.
      enddo
      topres=t
      return
      end
C
      subroutine readpdb(ipdb,opdb,ichain,natom,maxat,
     $     name,resn,chain,resid,x,y,z,wmain,newch,nchain,
     $     filnam,type,q19)
C
      implicit none
C     
C     Description of variables:
C==================================================================
C name(s)  max    actual Description
C          size   size   
C                        
C watn       1      1    Name of segmet for waters, starts with CWAT (C*4)
C wati       1      1    counter for water segments (I)
C line       1      1    a line from PDB or other file formats(C*132)
C
      integer ipdb,opdb,ichain,natom,maxat,nchain(*)
      real*8 x(*),y(*),z(*),wmain(*)
      character*(*) name(*),resn(*),chain(*),resid(*),newch(*)
      character*(*) type(*),filnam
      logical q19,dnaflag
C
      character*132 line
      character*100 tmp
      character*2 cc
      character*4 watn
C
      logical realatom
      integer i,lcc,lline,ion,terms,wati
      integer linlen
C
      wati=0
      watn(1:4)='CWAT'
C
      call letts(ichain,cc,lcc)
      write(tmp,'(a,''-'',a,''.PDB'')')filnam,cc(1:lcc)
      open(opdb,file=tmp,status='unknown')
C
C-----------------------------------
 100  read(ipdb,'(a)', end=101)line
      lline=linlen(line)
C
C     First make some IUPAC/CHARMM name adjustments.
C
      if (((line(1:4).eq.'ATOM').or.(line(1:6).eq.'HETATM'))
     $     .and.(.not.q19)) then
         if(line(18:20).eq.'HIS') line(18:20)='HSD'
         if(line(18:20).eq.'HOH') line(18:21)='TIP3'
         if((line(18:21).eq.'TIP3').and.(line(14:17).eq.'O   '))
     $        line(14:17)='OH2 '
         if((line(18:20).eq.'ILE').and.(line(14:17).eq.'CD1 '))
     $        line(14:17)='CD  '
C
C     The following is for RNA/DNA
C
         dnaflag=.false.
         if(line(18:20).eq.'  A') then
            line(18:20)='ADE'
            dnaflag=.true.
         endif
         if(line(18:20).eq.'  C') then
            line(18:20)='CYT'
            dnaflag=.true.
         endif
         if(line(18:20).eq.' +C') then
            line(18:20)='CYT'
            dnaflag=.true.
         endif
         if(line(18:20).eq.'  G') then
            line(18:20)='GUA'
            dnaflag=.true.
         endif
         if(line(18:20).eq.'  T') then
            line(18:20)='THY'
            dnaflag=.true.
         endif

         if(dnaflag.and.line(16:16).eq.'*') line(16:16)=' '

      endif
C
      realatom=.false.
      if(line(1:3).eq.'TER')terms=terms+1
      if(line(1:4).eq.'ATOM') then
         realatom=.true.
         natom=natom+1
         read(line,'(12x,a4,1x,a4,a1,1x,a4,3x,3f8.3,6x,f6.2)')
     $        name(natom),resn(natom),chain(natom),resid(natom),
     $        x(natom),y(natom),z(natom),wmain(natom)
         call strcln(name(natom))
         call namcln(name(natom))
         call strcln(resn(natom))
         call strcln(resid(natom))
         call settyp(1,name(natom),type(natom))
      endif
      if(line(1:6).eq.'HETATM') then
         realatom=.true.
         natom=natom+1
         read(line,'(12x,a4,1x,a4,a1,1x,a4,3x,3f8.3,6x,f6.2)')
     $        name(natom),resn(natom),chain(natom),resid(natom),
     $        x(natom),y(natom),z(natom),wmain(natom)
         call strcln(name(natom))
         call namcln(name(natom))
         call strcln(resn(natom))
         call strcln(resid(natom))
         call settyp(1,name(natom),type(natom))
         if(natom.gt.1)then
            if ((resn(natom).ne.resn(natom-1)).and.
     $           ((resn(natom).eq.'TIP3').or.(resn(natom).eq.'HOH ')))
     $           then
               wati=wati+1
               if(wati.eq.1)then
                  watn(1:4)='CWAT'
               else
                  call letts(wati,cc,lcc) 
                  write(watn,'(''CW'',a)')cc(1:lcc)
               endif
            endif
         endif
         if((resn(natom).eq.'TIP3').or.(resn(natom).eq.'HOH ')) then
            chain(natom)=watn(1:4)
            call settyp(0,name(natom),type(natom))
         endif
         if(chain(natom)(1:1).eq.' ') then
C
C     CHECK:
C           Defer the assignement of the IONs to the next
C           atom read. If only one atom in a residue then
C           assign as ION, else it is some other chain name:
C           name it with two or more caracters - so it doesn't
C           coincides with the other which may come later...
C
C     CHECK:
C           Do we want all the residues in empty (' ') chain as
C           one segment - if yes then make a flag for this!!
C
C
            call settyp(1,name(natom),type(natom))
C
            if(natom.gt.1) then
               if(resid(natom-1).ne.resid(natom)) then
                  if(ion.lt.10) then
                     write(chain(natom),'(''MOL'',i1)')ion
                     ion=ion+1
                  else if(ion.lt.100) then
                     write(chain(natom),'(''MO'',i2)')ion
                     ion=ion+1
                  else
                     write(*,*)'Fatal error: Too many molecules!!'
                     stop
                  endif
               else
                  chain(natom) = chain(natom-1)
               endif
            endif
         endif
      endif
C
C     Write each chain to a separate PDB file! This is because
C     we want to keep the resid names as originally in the PDB file!
C     In order to properly read PDB in CHARMM each chain/segment has
C     to be in a separate file. We also provide a sequence statements
C     (commented out) for further manipulation of the system by hand.
C
      if (realatom) then
         if (natom.gt.1) then
            if(chain(natom).ne.chain(natom-1)) then
               ichain=ichain+1
               write(opdb,'(''TER'')')
               write(opdb,'(''END'')')
               call letts(ichain,cc,lcc)
               write(tmp,'(a,''-'',a,''.PDB'')')
     $              filnam,cc(1:lcc)
               open(opdb,file=tmp,status='unknown')
               write(opdb,'(a)')line(1:lline)
               nchain(ichain)=natom
            else
               write(opdb,'(a)')line(1:lline)
            endif
         else
            write(opdb,'(a)')line(1:lline)
         endif
C
C        Check here if the chain name was already used before
C        Some PDB files have more than 36 chains, so they start
C        to repeat. In this case add another letter to the PDB
C        chain name. This has to be done for every atom!
C
CC         write(*,*)'natom,ichain,nchain,chain=',
CC     $        natom,ichain,nchain(ichain),chain(natom)
         call chkchain(natom,ichain,nchain,chain,newch)
C
C
      endif
C
      if (natom.ge.maxat) then
         write(*,*)'Too many atoms: modify the program!'
         stop
      endif
C
      goto 100
C
 101  continue
C
C     TODO: This is not very well thought out. Try to fix!
      do i=1,natom
         chain(i)=newch(i)
      enddo
C
C     End of main READ section
C
      return
      end
C
      subroutine readcrd(icrd,opdb,oxyz,filnam,ichain,
     $     x,y,z,wmain,title,ntitle,natom,maxat,
     $     resn,name,chain,resid,type,nchain)
C
      implicit none
C
C     This is CHARMM native coordinate format reader
C
C     Here we assume everything is pretty much OK
C     But we need to distinguish between residues for which we already
C     have topology, and the rest (there is no HETATM command in .crd
C     files!)  We also need to preserve resid names, which means we
C     still have to go through PDB files for each chain!!!. As a bonus
C     we have then a .crd to PDB converter!!! The one in CHARMM is not
C     that good - it ignores the chain column among other things...
C
C
C     Description of variables:
C==================================================================
C name(s)  max    actual Description
C          size   size   
C                        
C lline      1      1    length of line
C line       1      1    a line from PDB or other file formats(C*132)
C
      integer icrd,opdb,oxyz,ichain,ntitle,natom,maxat,nchain(*)
      real*8 x(*),y(*),z(*),wmain(*)
      character*(*) filnam,title(*),resn(*),name(*),chain(*)
      character*(*) type(*),resid(*)
C
      integer i,j,lline,linlen,ires
      character*100 tmp
      character*132 line
      logical topres
      character*6 pdbcmd
      character*4 pdbato,pdbrid
      character*1 artchn(26),tmpchn
      real*8 one
C
      data one/1.0d0/
C
C     This is to automatically assign chain names for PDB
C     TODO: Better use `call letts'
      do i = 1, 26
         artchn(i)=char(ichar('A')+i-1)
      enddo
C
      write(tmp,'(a,''-'',i1,''.PDB'')')filnam,ichain
      open(opdb,file=tmp,status='unknown')
C
C-----------------------------------
C
C     First read the title and save for output
C
      i=1
 200  read(icrd,'(a)', end=201)line
      lline=linlen(line)
      if (line(1:1).eq.'*') then
         title(i)=line(1:lline)
         i=i+1
         goto 200
      else
         goto 201
      endif
 201  continue
      ntitle=i-1
      read(line,'(i5)')natom
C
      if (natom.ge.maxat) then
         write(*,*)'Too many atoms: modify the program!'
         stop
      endif
C
      do i = 1, natom
         read(icrd,'(i5,i5,1x,a4,1x,a4,3f10.5,1x,a4,1x,a4,f10.5)')
     $        j,ires, resn(i), name(i), x(i), y(i), z(i), chain(i),
     $        resid(i),wmain(i)
         call strcln(name(i))
         call strcln(resn(i))
         call strcln(resid(i))
         if (topres(resn(i))) then
            pdbcmd='ATOM  '
            if(resn(i).eq.'TIP3')pdbcmd='HETATM'
            if(resn(i).eq.'HOH ')pdbcmd='HETATM'
            call settyp(0,name(i),type(i))
         else
            pdbcmd='HETATM'
            call settyp(1,name(i),type(i))
         endif
C
C     Write each chain to a separate PDB file! This is because
C     we want to keep the residue IDs as originally in the PDB file!
C     In order to properly read PDB in CHARMM each chain/segment has
C     to be in a separate file. We also provide a sequence statements
C     (commented out) for further manipulation of the system by hand.
C
         if (i.gt.1) then
            if(chain(i).ne.chain(i-1)) then
               ichain=ichain+1
               write(opdb,'(''TER'')')
               write(opdb,'(''END'')')
               write(tmp,'(a,''-'',i1,''.PDB'')')
     $              filnam,ichain
               open(opdb,file=tmp,status='unknown')
               nchain(ichain)=i
            endif
         endif
         tmpchn=chain(i)(1:1)
         if (linlen(chain(i)).gt.1) tmpchn=artchn(ichain)
         pdbato=name(i)
         if(linlen(name(i)).lt.4)pdbato=' '//name(i)(1:3)
         pdbrid=resid(i)
         call strlj(pdbrid)
         write(opdb,
     $        '(a6,i5,1x,a4,1x,a4,a1,1x,a4,3x,3f8.3,2f6.2,6x,a4,1x,a1)')
     $        pdbcmd,i,pdbato,resn(i),tmpchn,pdbrid,
     $        x(i),y(i),z(i),one,wmain(i),chain(i),name(i)(1:1)
         write(oxyz,
     $        '(a6,i5,1x,a4,1x,a4,a1,1x,a4,3x,3f8.3,2f6.2,6x,a4,1x,a1)')
     $        pdbcmd,i,pdbato,resn(i),tmpchn,pdbrid,
     $        x(i),y(i),z(i),one,wmain(i),chain(i),name(i)(1:1)
C
      enddo
C
      return
      end
C
      logical function cblack(c)
C
C     find out if ths character c is white space or delimiter
C     or something else
C
      implicit none
      character*1 c
C
      cblack=.true.
      if(c.eq.' ')cblack=.false.
      if(ichar(c).eq.9)cblack=.false.
C
      return
      end
C
      subroutine upcase(s)
C
      implicit none
      character*(*) s
      integer k,i,l,linlen,d
C
      l=linlen(s)
      d=ichar('a')-ichar('A')
      do i=1,l
         k=ichar(s(i:i))
         if(k.ge.ichar('a')) s(i:i)=char(k-d)
      enddo
C
      return
      end
C
      integer function nword(l)
C
      implicit none
      character*(*) l
      integer i,k,f,linlen
      logical cblack
C
      k=linlen(l)
      nword=0
      if(k.le.0) return
      f=1
      do i=1,k
         if((f.eq.1).and.cblack(l(i:i)))then
            nword=nword+1
            f=0
         else if(l(i:i).eq.' ') then
            f=1
         endif
      enddo
C
      return
      end
C
      subroutine readxyz(ixyz,ocrd,opdb,natom,maxat,filnam,
     $     x,y,z,wmain,symb,name,resn,resid,chain,type,maxyzres)
C
      implicit none
C
C+++++++++++++++++++++++++++++++++++++
C
C     Read XYZ file: Here we assume that everything goes into one
C     segment, and we make artificial residues
C     (some number of atoms = maxyzres)
C     in order not to upset CHARMM too much.
C     residue names = first 2 characters from the file name
C     segment name = first 4 characters from the file name
C     atom name = XYZ accepts only periodic table names, so we add
C                 numbers and characters to it
C
C     
C     Description of variables:
C==================================================================
C name(s)  max    actual Description
C          size   size   
C
C line       1      1    a line from PDB or other file formats(C*132)
C                        
      integer ixyz,ocrd,opdb,natom,maxat,maxyzres
      real*8 x(*),y(*),z(*),wmain(*)
      character*(*) filnam,symb(*),name(*),resn(*),resid(*)
      character*(*) chain(*),type(*)
C
      integer lcc,i,j,k,l,ires,linlen,nword
      real*8 one,zero
      character*132 line,tmp
      character*4 pdbrid,pdbato
      character*2 cc
      character*6 pdbcmd
C
      data one/1.0d0/,zero/0.0d0/
C
      write(tmp,'(a,''-A.PDB'')')filnam
      open(15,file=tmp,status='unknown')
C
C     Start of the main READ block:
C
C-----------------------------------
C
C     First read the title and save for output
C
      read(ixyz,*)natom
      read(ixyz,'(a)')line
C
      if (natom.ge.maxat) then
         write(*,*)'Too many atoms: modify the program!'
         stop
      endif
C
      lcc=linlen(filnam)
      write(ocrd,
     $     '(''* generated coordinate file from XYZ format: '',a)')
     $     filnam(1:lcc)
      lcc=linlen(line)
      if(lcc.gt.0)write(ocrd,'(''* '',a)')line(1:lcc)
      write(ocrd,'(''*'')')
C
C     The rest of coordinate (.crd) file is written elsewhere
C     
      do i = 1, natom
         read(ixyz,'(a)')line
         l=nword(line)
         if(l.ge.5)then
            read(line,*)symb(i),x(i),y(i),z(i),wmain(i)
         else if(l.eq.4) then
            wmain(i)=zero
            read(line,*)symb(i),x(i),y(i),z(i)
         else
            write(*,*)'This file may not be the normal xyz file'
         endif
C
         call strcln(symb(i))
         pdbcmd='HETATM'
         call settyp(1,symb(i),type(i))
C
C     Construct the necessary information for CRD and PDB formats
C     Chain name is always first letters (upto 4) from the filename
         lcc=linlen(filnam)
         chain(i)='    '
         chain(i)=filnam(1:min(4,lcc))
         call upcase(chain(i))
C     Atom name is atomic symbol, plus number, plus upto 2 letter code
C     To make this better we could construct the atoms like:
C     First only the symbol name. If this is already taken then add one number
C     So counting for each symbolic name can be separate. For now we just
C     do the simpler version. Always add a number and if necessary also 
C     letters from A - ZZ. Don't distinguish between different symbols
C     so this code is good only for about 6000 atoms
C
C     FIXME: This code doesn't work very well !! ???
C
         j=linlen(symb(i))
         if(i.lt.10) then
            write(pdbato,'(a,i1)')symb(i)(1:j),i
         else
            k=(i-10)/26+1
            if(k.le.9) then
               l=i-9-(k-1)*26
               call letts(l,cc,lcc)
               write(pdbato,'(a,i1,a)')symb(i)(1:j),k,cc(1:lcc)
            else
               k=(i-10-26*9)/(26*26)+1
               l=i-9-26*8-(k-1)*26*26
               call letts(l,cc,lcc)
               write(pdbato,'(a,i1,a)')symb(i)(1:j),k,cc(1:lcc)
            endif
         endif
         name(i)=pdbato
         ires=i/maxyzres+1
         call letts(ires,cc,lcc)
         resn(i)=cc(1:lcc)
C
C     For residue names we want something better than just ABC...
C
         lcc=linlen(filnam)
         resn(i)='    '
         resn(i)=filnam(1:min(4,lcc))
         call upcase(resn(i))
C
         write(resid(i),'(i4)')ires
C
         if(linlen(name(i)).lt.4)pdbato=' '//name(i)(1:3)
         pdbrid=resid(i)
         call strlj(pdbrid)
         write(opdb,
     $        '(a6,i5,1x,a4,1x,a4,a1,1x,a4,3x,3f8.3,2f6.2,6x,a4,1x,a1)')
     $        pdbcmd,i,pdbato,resn(i),chain(i),pdbrid,
     $        x(i),y(i),z(i),one,wmain(i),chain(i),name(i)(1:1)
C
      enddo
C
      return
      end
C
      subroutine readmol(imol,ocrd,opdb,natom,nbonds,maxat,filnam,
     $     x,y,z,wmain,symb,ibond,jbond,name,resn,resid,chain,
     $     type,maxyzres)
C
      implicit none
C
C+++++++++++++++++++++++++++++++++++++
C
C     Read MOL file (not MOL2):
C     Here we assume that everything goes
C     into one segment, and we make artificial residues
C     (some number of atoms = maxyzres)
C     in order not to upset CHARMM too much.
C     residue names = first 2 characters from the file name
C     segment name = first 4 characters from the file name
C     atom name = MOL accepts only periodic table names, so we add
C                 numbers and characters to it
C
C     
C     Description of variables:
C==================================================================
C name(s)  max    actual Description
C          size   size   
C
C line       1      1    a line from PDB or other file formats(C*132)
C                        
      integer imol,ocrd,opdb,natom,nbonds,maxat,maxyzres
      real*8 x(*),y(*),z(*),wmain(*)
      integer ibond(*),jbond(*)
      character*(*) filnam,symb(*),name(*),resn(*),resid(*)
      character*(*) chain(*),type(*)
C
      integer lcc,i,j,k,l,ires,linlen,nword
      real*8 one,zero
      character*132 line,tmp
      character*4 pdbrid,pdbato
      character*2 cc
      character*6 pdbcmd
C
      data one/1.0d0/,zero/0.0d0/
C
      write(tmp,'(a,''-A.PDB'')')filnam
      open(15,file=tmp,status='unknown')
C
C     Start of the main READ block:
C
C-----------------------------------
C
C     First read the title and save for output
C
      read(imol,'(a)')line
      read(imol,*)
      read(imol,*)
      read(imol,'(2i3)')natom,nbonds
C
      if (natom.ge.maxat) then
         write(*,*)'Too many atoms: modify the program!'
         stop
      endif
C
      lcc=linlen(filnam)
      write(ocrd,
     $     '(''* generated coordinate file from MOL format: '',a)')
     $     filnam(1:lcc)
      lcc=linlen(line)
      if(lcc.gt.0)write(ocrd,'(''* '',a)')line(1:lcc)
      write(ocrd,'(''*'')')
C
C     The rest of coordinate (.crd) file is written elsewhere
C     
      do i = 1, natom
         read(imol,'(3f10.4,1x,a2)')x(i),y(i),z(i),symb(i)
         wmain(i)=zero
C
         call strcln(symb(i))
         pdbcmd='HETATM'
         call settyp(1,symb(i),type(i))
C
C     Construct the necessary information for CRD and PDB formats
C     Chain name is always first letters (upto 4) from the filename
         lcc=linlen(filnam)
         chain(i)='    '
         chain(i)=filnam(1:min(4,lcc))
         call upcase(chain(i))
C     Atom name is atomic symbol, plus number, plus upto 2 letter code
C     To make this better we could construct the atoms like:
C     First only the symbol name. If this is already taken then add one number
C     So counting for each symbolic name can be separate. For now we just
C     do the simpler version. Always add a number and if necessary also 
C     letters from A - ZZ. Don't distinguish between different symbols
C     so this code is good only for about 6000 atoms
         j=linlen(symb(i))
         if(i.lt.10) then
            write(pdbato,'(a,i1)')symb(i)(1:j),i
         else
            k=(i-10)/26+1
            if(k.le.9) then
               l=i-9-(k-1)*26
               call letts(l,cc,lcc)
               write(pdbato,'(a,i1,a)')symb(i)(1:j),k,cc(1:lcc)
            else
               k=i/676+1
               l=i-10-(k-1)*676
               call letts(l,cc,lcc)
               write(pdbato,'(a,i1,a)')symb(i)(1:j),k,cc(1:lcc)
            endif
         endif
         name(i)=pdbato
         ires=i/maxyzres+1
         call letts(ires,cc,lcc)
         resn(i)=cc(1:lcc)
C
C     For residue names we want something better than just ABC...
C
         lcc=linlen(filnam)
         resn(i)='    '
         resn(i)=filnam(1:min(4,lcc))
         call upcase(resn(i))
C
         write(resid(i),'(i4)')ires
C
         if(linlen(name(i)).lt.4)pdbato=' '//name(i)(1:3)
         pdbrid=resid(i)
         call strlj(pdbrid)
         write(opdb,
     $        '(a6,i5,1x,a4,1x,a4,a1,1x,a4,3x,3f8.3,2f6.2,6x,a4,1x,a1)')
     $        pdbcmd,i,pdbato,resn(i),chain(i),pdbrid,
     $        x(i),y(i),z(i),one,wmain(i),chain(i),name(i)(1:1)
C
      enddo
C
C     MOL file has also info about the bonds...
C
      if(nbonds.gt.maxat)then
         write(*,*)'MOL reader: too many bonds!'
         return
      endif
      do k = 1, nbonds
         read(imol,'(2i3)')i,j
         ibond(k)=i
         jbond(k)=j
      enddo
C
      return
      end
C
      subroutine readspdb(imol,ocrd,opdb,natom,nbonds,maxat,filnam,
     $     x,y,z,wmain,symb,ibond,jbond,name,resn,resid,chain,
     $     type,maxyzres,qcon)
C
      implicit none
C
C+++++++++++++++++++++++++++++++++++++
C
C     Read XYZ file: Here we assume that everything goes into one
C     segment, and we make artificial residues
C     (some number of atoms = maxyzres)
C     in order not to upset CHARMM too much.
C     residue names = first 2 characters from the file name
C     segment name = first 4 characters from the file name
C     atom name = XYZ accepts only periodic table names, so we add
C                 numbers and characters to it
C
C     
C     Description of variables:
C==================================================================
C name(s)  max    actual Description
C          size   size   
C
C line       1      1    a line from PDB or other file formats(C*132)
C                        
      integer imol,ocrd,opdb,natom,nbonds,maxat,maxyzres
      real*8 x(*),y(*),z(*),wmain(*)
      integer ibond(*),jbond(*)
      character*(*) filnam,symb(*),name(*),resn(*),resid(*)
      character*(*) chain(*),type(*)
      logical qcon
C
      integer lcc,i,j,k,l,ires,linlen,nword
      real*8 one,zero
      character*132 line,tmp
      character*4 pdbrid,pdbato
      character*2 cc
      character*6 pdbcmd
C
      data one/1.0d0/,zero/0.0d0/
C
      write(tmp,'(a,''-A.PDB'')')filnam
      open(15,file=tmp,status='unknown')
C
C     Start of the main READ block:
C
C-----------------------------------
C
C     First read the title and save for output
C
      read(imol,'(a)')line
      read(imol,*)
      read(imol,*)
      read(imol,'(2i3)')natom,nbonds
C
      if (natom.ge.maxat) then
         write(*,*)'Too many atoms: modify the program!'
         stop
      endif
C
      lcc=linlen(filnam)
      write(ocrd,
     $     '(''* generated coordinate file from XYZ format: '',a)')
     $     filnam(1:lcc)
      lcc=linlen(line)
      if(lcc.gt.0)write(ocrd,'(''* '',a)')line(1:lcc)
      write(ocrd,'(''*'')')
C
C     The rest of coordinate (.crd) file is written elsewhere
C     
      do i = 1, natom
         read(imol,'(3f10.4,1x,a2)')x(i),y(i),z(i),symb(i)
         wmain(i)=zero
C
         call strcln(symb(i))
         pdbcmd='HETATM'
         call settyp(1,symb(i),type(i))
C
C     Construct the necessary information for CRD and PDB formats
C     Chain name is always first letters (upto 4) from the filename
         lcc=linlen(filnam)
         chain(i)='    '
         chain(i)=filnam(1:min(4,lcc))
         call upcase(chain(i))
C     Atom name is atomic symbol, plus number, plus upto 2 letter code
C     To make this better we could construct the atoms like:
C     First only the symbol name. If this is already taken then add one number
C     So counting for each symbolic name can be separate. For now we just
C     do the simpler version. Always add a number and if necessary also 
C     letters from A - ZZ. Don't distinguish between different symbols
C     so this code is good only for about 6000 atoms
         j=linlen(symb(i))
         if(i.lt.10) then
            write(pdbato,'(a,i1)')symb(i)(1:j),i
         else
            k=(i-10)/26+1
            if(k.le.9) then
               l=i-9-(k-1)*26
               call letts(l,cc,lcc)
               write(pdbato,'(a,i1,a)')symb(i)(1:j),k,cc(1:lcc)
            else
               k=i/676+1
               l=i-10-(k-1)*676
               call letts(l,cc,lcc)
               write(pdbato,'(a,i1,a)')symb(i)(1:j),k,cc(1:lcc)
            endif
         endif
         name(i)=pdbato
         ires=i/maxyzres+1
         call letts(ires,cc,lcc)
         resn(i)=cc(1:lcc)
C
C     For residue names we want something better than just ABC...
C
         lcc=linlen(filnam)
         resn(i)='    '
         resn(i)=filnam(1:min(4,lcc))
         call upcase(resn(i))
C
         write(resid(i),'(i4)')ires
C
         if(linlen(name(i)).lt.4)pdbato=' '//name(i)(1:3)
         pdbrid=resid(i)
         call strlj(pdbrid)
         write(opdb,
     $        '(a6,i5,1x,a4,1x,a4,a1,1x,a4,3x,3f8.3,2f6.2,6x,a4,1x,a1)')
     $        pdbcmd,i,pdbato,resn(i),chain(i),pdbrid,
     $        x(i),y(i),z(i),one,wmain(i),chain(i),name(i)(1:1)
C
      enddo
C
C     MOL file has also info about the bonds...
C
      if(nbonds.gt.maxat)then
         write(*,*)'MOL reader: too many bonds!'
         return
      endif
      do k = 1, nbonds
         read(imol,'(2i3)')i,j
         ibond(k)=i
         jbond(k)=j
      enddo
C
      return
      end
C
      subroutine testlet
      integer i,l
      character*5 c
C
      open(unit=87,file='letts',status='unknown')
      do i=1,3000
         call letts(i,c,l)
         write(87,*)'i,l=',i,l,':',c(1:l)
      enddo
      return
      end
